beagle<-read.table("beagle_r2.gz",sep="\t")
impute<-read.table("impute_r2.gz",sep="\t")
colnames(beagle)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
colnames(impute)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
dim(beagle)
dim(impute)
pos<-read.table("position",sep=" ")
dim(pos)
pos[1,]
pos<-read.table("position",sep=" ",head=T)
pos[1,]
merged<-merge(beagle,impute,by="marker_identifier")
colnames(pos)<-c("position","marker_identifier")
merged<-merge(merged,pos,by="marker_identifier")
dim(merged)
merged[1,]
dim(beagle)
dim(impute)
pdf("snps_dos_beagle_vs_impute.pdf")
plot(merged$dosage_r_squared.x,merged$dosage_r_squared.y,xlab = c("beagle dosage r2"),ylab=c("impute dosage r2"))
cor(merged$dosage_r_squared.x,merged$dosage_r_squared.y,method="pearson")
main<-paste("Beagle dos r2 VS Impute dos r2 ,pearson corr=",cor(merged$dosage_r_squared.x,merged$dosage_r_squared.y,method="pearson"),sep=" ")
dev.off()
pdf("snps_dos_beagle_vs_impute.pdf")
plot(merged$dosage_r_squared.x,merged$dosage_r_squared.y,xlab = c("beagle dosage r2"),ylab=c("impute dosage r2"))
dev.off()
pdf("snps_dos_beagle_vs_impute.pdf")
plot(merged$dosage_r_squared.x,merged$dosage_r_squared.y,xlab = c("beagle dosage r2"),ylab=c("impute dosage r2"),pch=".",col=c("red"))
main<-paste("Beagle dos r2 VS Impute dos r2 ,pearson corr=",cor(merged$dosage_r_squared.x,merged$dosage_r_squared.y,method="pearson"),sep=" ")
dev.off()
pdf("snps_dos_beagle_vs_impute.pdf")
main<-paste("Beagle dos r2 VS Impute dos r2 ,pearson corr=",cor(merged$dosage_r_squared.x,merged$dosage_r_squared.y,method="pearson"),sep=" ")
main
plot(merged$dosage_r_squared.x,merged$dosage_r_squared.y,xlab = c("beagle dosage r2"),ylab=c("impute dosage r2"),pch=".",col=c("red"))
abline(lm(merged$dosage_r_squared.y~merged$dosage_r_squared.x))
dev.off()
pdf("snps_dos_beagle_vs_impute.pdf")
main<-paste("Beagle DosR2 vs ImputeDosR2 ,pearson corr=",cor(merged$dosage_r_squared.x,merged$dosage_r_squared.y,method="pearson"),sep=" ")
main
 plot(merged$dosage_r_squared.x,merged$dosage_r_squared.y,xlab = c("beagle dosage r2"),ylab=c("impute dosage r2"),pch=".",col=c("red"),main=main)
abline(lm(merged$dosage_r_squared.y~merged$dosage_r_squared.x))
dev.off()
pdf("Position.pdf")
merged[1,]
plot(merged$position,merged$dosage_r_squared.x,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
points(merged$position,merged$dosage_r_squared.y,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("green"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
dev.off()
pdf("Position.pdf",width=20,height=10)
plot(merged$position,merged$dosage_r_squared.x,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
points(merged$position,merged$dosage_r_squared.y,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("green"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
dev.off()
pdf("Position.pdf",width=200,height=10)
plot(merged$position,merged$dosage_r_squared.x,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
points(merged$position,merged$dosage_r_squared.y,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("green"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
dev.off()
pdf("Position.pdf",width=500)
plot(merged$position,merged$dosage_r_squared.x,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
points(merged$position,merged$dosage_r_squared.y,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("green"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
dev.off()
pdf("Position.pdf",width=200)
plot(merged$position,merged$dosage_r_squared.x,xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
points(merged$position,merged$dosage_r_squared.y,pch=".",col=c("green"))
dev.off()

pdf("Position.pdf",width=200)
start<-c(16894264,25481522,34068780,42656038)
stop<-c(25481522,34068780,42656038,51243298)
i<-1
while(i<=length(start))
{
	plot(merged[which(merged[,16]>=start[i] & merged[,16]<stop[i]),"position"],merged[which(merged[,16]>=start[i] & merged[,16]<stop[i]),"dosage_r_squared.x"],xlab = c("POSITION"),ylab=c("Dosage r2"),pch=".",col=c("red"),main="Beagle Dos R2(red) vs Impute Dos R2(green)")
	points(merged[which(merged[,16]>=start[i] & merged[,16]<stop[i]),"position"],merged[which(merged[,16]>=start[i] & merged[,16]<stop[i]),"dosage_r_squared.y"],pch=".",col=c("green"))
	i<-i+1
}
dev.off()
multhist <- function(x,main,x.loc,...) {
  new.x <- lapply(x,hist,plot=F)
  tmp <- lapply(new.x,function(x) x$breaks)
  max.break <- lapply(tmp,function(x) max(x))
  max.break <- which(unlist(max.break)==max(unlist(max.break)))[1]
  breaks <- new.x[[max.break]]$breaks
  new.x <- lapply(x,hist,plot=F,breaks=breaks)
  y.max <- max(unlist(lapply(new.x,function(x) max(x$counts))))
  barplot(do.call(rbind,lapply(new.x,function(x) x$counts)),beside=T,
          space=c(0,0),names.arg=signif(new.x[[1]]$mids,2),
          axes=T,legend.text=T,main=main,args.legend=list(x=x.loc,y=y.max),...)
}
pdf("Multiple_hist_chr22.pdf")
merged[1,]
merged[1:2,5]
merged[1:2,5:10]
merged[1:2,5:12]
multhist(x=list(Beagle_r2= merged_f[,5],Impute_r2= merged_f[,12]),main="Beagle_r2 VS Impute_r2", x.loc=20)
multhist(x=list(Beagle_r2= merged[,5],Impute_r2= merged[,12]),main="Beagle_r2 VS Impute_r2", x.loc=20)
dev.off()
merged[1,]
max(merged[,1])
dim(merged)
max(merged[,16])
min(merged[,16])
getwd()
51243297-16894264
(51243297-16894264)/4
16894264+8587258
25481522+8587258
34068780+8587258
42656038+8587258
getwd()
savehistory("newplots_script.R")
